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We consider classical hard-core particles moving on two parallel chains in the same direction. An 
interaction between the channels is included via the hopping rates. For a ring, the stationary state 
has a product form. For the case of coupling to two reservoirs, it is investigated analytically and 
numerically. In addition to the known one-channel phases, two new regions are found, in particular 
one, where the total density is fixed, but the filling of the individual chains changes back and 
forth, with a preference for strongly different densities. The corresponding probability distribution 
is determined and shown to have a universal form. The phase diagram and general aspects of the 
problem are discussed. 



I. INTRODUCTION 

Driven many-particle systems have been the topic of numerous studies in recent years |^ . The simplest example is a 
one-dimensional lattice gas where hard-core particles move stochastically in one direction. This model, also known as 
asymmetric exclusion process, can be treated exactly and therefore has become a reference system in this field 
The characteristic feature of such driven systems is a nonzero current in the stationary state. For the case of open 
boundaries, this current transports the 'information' from the boundaries into the bulk and nonequilibrium phase 
transitions can arise [^,^. These do not have analogues in equilibrium systems where the boundaries normally do 
not play a significant role. It was shown in |0] that the phase transitions for a generic driven system with one type of 
particles are governed entirely by an extremal principle for the macroscopic current j{p), where p denotes the average 
density of the particles. This principle states that the stationary bulk current assumes its absolute minimum within 
the interval set by the boundary densities, if /9_(left) < p+ (right) and the flux is towards the right. If /3_ > it 
takes its absolute maximum. In either case, the current and the density fluctuate only slightly around their stationary 
values and the fluctuations vanish as the system size goes to infinity. 

In the present communication we study a system with two parallel channels and show that these two properties 
do not hold in general. Our example is a simple extension of the one-channel model and contains one additional 
parameter which measures the coupling between the channels. There is no exchange of particles, but the hopping 
rates in one chain depend on the local configuration in the other one. They are chosen in such a way that, for a ring, 
the stationary state has a simple product form and thus the current can be obtained explicitly. With reservoirs at the 
ends, the system was studied by a combination of numerical and analytical methods. In a - phase diagram, it 
shows a number of regions with equal densities in both chains and various values for the current. In addition, however, 
there are two other regions with unexpected new features. They develop out of the first-order transition line of the 
single-chain problem as the interaction between the channels is turned on. In both of them, the overall densities pi 
and p2 in the two chains fluctuate strongly and only a probability distribution w{pi,p2) can be given. Together with 
that, symmetry-breaking phenomena appear. 

In one region, one finds spatial coexistence between sections of equal and of unequal (but fixed) densities in the 
two chains, with the size of these sections varying in time. This is similar to the situation on the transition line for 
one chain, where sections of high and low density, separated by a domain wall, coexist. It is related to the fact that 
the current, which must be constant throughout the system (and thus plays the role of a chemical potential), can be 
the same for different densities. This will be called the mixed-phase region in the following. 

In the other region, which appears if the coupling exceeds a critical value and then grows at the expense of the 
first one, only the unsymmetric configurations exist. The total system is then practically half filled, pi + P2 = 1, 
but the individual densities change in time. The same holds for the current. The most probable configurations are 
those where one channel is relatively empty and the other relatively full. Between them the system diffuses back and 
forth and hence we will use the term seesaw region in the following. The time in which the channels interchange 
roles increases only as a power of the system size, in contrast to the result for a model with two kinds of particles on 
one chain 1^,^ or simplified versions of it ||9|,[l0| . The situation can be compared to that at an equilibrium first-order 
transition with a vanishing or size-independent free-energy barrier between the phases. In this sense, the symmetry 
breaking could be called 'weak'. The origin of this behaviour is related to the existence of fast and slow processes in 
the system, as will be discussed in detail. 
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The paper is organized as follows. In section II we define the model and describe the solution on a ring. In section 
III we treat it with boundary reservoirs and present numerical results from Monte-Carlo and mean-field calculations. 
The case of strong interactions, where one can simplify the problem and obtain analytical results for the seesaw region, 
is the topic of section IV. After that we turn in section V to general interactions and establish the complete — p-t^ 
phase diagram. Section VI contains the conclusion and some discussion of open problems. Some details concerning 
the boundary rates and the mean-field equations can be found in Appendices A and B. 



II. MODEL 



Our model consists of two parallel chains, on each of which particles can hop towards the right if the next site is 
empty. The hopping rate in one chain depends on the configuration of the neighbouring sites in the other and one 
has the four processes shown in Fig. |l|. The stationary state has a simple form (given below) if the rates satisfy the 
condition a -\- (3 = 27. In the following we choose a — (3 — ^— \,so that the remaining rate 

e = exp{—v) (1) 

is the only parameter. It refers to those processes where the site besides the jumping particle is empty, but the next 
one in the forward direction is occupied. We will always assume e < 1. This can be viewed as the result of short-range 
interactions which increase the barrier similar as in pT| . In a traffic- model context, e would describe a hesitation to 
move besides another car. 




FIG. 1. The four elementary hopping processes, shown here for the first chain, and their rates. In the study, the first three 
rates are set equal to 1. 



The state of the system can be described either by occupation numbers nk,mk or by spin variables crk,Tk for the 
two chains. The first take the values 1(0), the latter the values -l-l(-l) if site k is occupied (empty). For a ring with 
iV sites, the stationary probability P{a, r) = P{(Ji, a2, ■ ■ ■ , crjv, n, T2, . . . , tn) then has the Boltzmann form 

^ 1 

P{<j,T)=Cl[exp{--,^<Jkn) (2) 
fc=i 

which means that different sites k are uncorrelated. This can be proved by a straightforward consideration of the gain 
and loss processes. For the ring geometry, the stationary densities are constant and the quantities of main interest 
are the currents ji, ^2 in the two chains. These can be calculated either by working in a grand canonical ensemble or 
via the mean-field equations of Appendix B which are exact here due to the form of (H). For example, the mean-field 



expression for ji is, from (Bl), with pk — nkirik 

ji = Ukil - Uk+i) + (e - l)(n/c -~Pk)imk+i ~ Pk+i) (3) 

where all quantities are expectation values. In the stationary state they are independent of k, Uk = pi,mk = P2 and 
one only has to determine p. The expression for the current for arbitrary densities pi , p2 is complicated and is omitted 
here. In particular cases pi — P2 pi = 1 — p2, the currents in the two chains are the same and given by 



1 ± [1 - ^TTf{±^)) /F{±v,p) 



2 



(4) 



F{v,p) = A{e- -\)p{\~ p), (5) 

where the upper (lower) sign holds for pi = P2 — P and pi — 1 — p2 — p, respectively. The resulting curves are shown 
in Fig. ^ for various values of the interaction parameter v. One can see that for larger v a minimum at p = 1/2 
exists in the case of equal densities. The value where this first happens is i^crit = InA w 1.39. The reason is that for 
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half filling and large v adjacent sites are mainly occupied by particle-hole pairs and then only hopping with the small 
rate e is possible. Such a double-peak structure of the current leads to a much richer phase diagram in the one-chain 
problem with reservoirs and will be important here, too. 



(a) (b) 
0.2 I ^ 1 0.2 I ^ 




FIG. 2. Current j vs. density p for a ring according to (3,4) for four values of the interaction, v = 1.2,2,4,8 (a-d). Upper 
and lower curves correspond to p = p\ = p2 and p = pi = 1 — /92, respectively. 

With reservoirs, particles enter the system at site n — 1 and leave from site N. If one adds reservoir sites n = 
and n = iV-l- 1 to the chains, the bulk processes of Fig. |l| appear also at the boundaries. The rates can then be chosen 
in such a way that the dynamics is the same as in the interior, but at prescribed densities p- and p+ These 
are thus the only boundary parameters which enter. The procedure is described in more detail in Appendix A. In 
the following we will normally choose equal boundary densities for both chains which makes the problem completely 
symmetric between them. Since the transport can be viewed as the motion of vacancies towards the left, it is also 
symmetric under the exchange /9_ ^ 1 — p+ , Pq, ^ 1 — Pa- 

III. NUMERICS 

We have performed Monte-Carlo (MC) simulations of our system for different values of p„, p+ and ly, looking first 
at global quantities, namely the overall densities pi{t) and p2{t) in the two chains. The motivation for this was that, 
for a single chain, the average density of particles is an order parameter characterizing the different phases. This 
means that p(t) fluctuates only slightly around its mean value and the fluctuations vanish as the system size grows. 

Typical results of such MC calculations are presented in Figj^ for 1^ — 2 and three different values of the boundary 
densities. While in Fig||(a) the behaviour is as for a single chain, the other two figures show a 'random walk' of pi 
and p2 within a large range of densities. These large fluctuations are intrinsic, they do not change qualitatively with 
the system size. In FigH(c) the situation is still relatively simple since the sum of the densities stays approximately 
constant. In Fig.||(b), however, the behaviour is rather irregular and the global densities do not give a sufficient 
description. Rather one has to look at the state of the system in more detail. It then turns out that the seemingly 
chaotic pattern is connected with the spatial coexistence of several phases. 
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FIG. 3. Time evolution of the average densities in the two chains for three typical cases, as obtained from simulations for 
a system with v = 2 and N = 1000 sites. Boundary densities P-,p+ : (a) 0.15, 0.5 ; (b) 0.15, 0.8 ; (c) 0.42, 0.6. 



This can already be seen in a mean-field (MF) analysis of the system. The corresponding equations are given in 
Appendix B. We integrated them over time, starting from random initial conditions and stopping when the current 
had converged to 10~^ in relative units. In the region of small p_ and large p+, the mean-field density profiles then 
had the typical shapes seen in FigJ|. 
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FIG. 4. Density profiles in the mean-field approximation for a system with v = 2 and A'^ = 200 sites. Boundary densities 
P-,p+ : (a) 0.15, 0.82; (b) 0.15, 0.85. The horizontal lines indicate the densities p,l — p predicted by equating the currents 
j{p,l — p) = j{p- , P-), see section III 
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Let us first discuss Fig.^(a) wiiicii corresponds to p_ < 1 — p+ . Here one lias two regions wliere the local densities 
are constant. On the left they are the same in both chains, {rik), (mk) = p-, while on the right they are different. 
This leads to the bubble-like structure in the figure. That such a coexistence is possible, follows already from Fig||, 
since a given value of the current (for the ring) can be realized in different ways. Correspondingly, the densities in the 
bubble region are given by (nk) — p, (mk) — 1 — p where j{p, 1 — p) — j{p-,p-). The left end of the bubble can be 
anywhere, its location depends on the initial conditions. Similarly, the two channels can exchange roles. However, the 
bubble is always 'glued' to the right boundary. For p_ > 1 — p+ one has a similar picture but the region with equal 
densities now appears on the right and the values are (rife), {mk) = p+. This is a consequence of the particle-hole 
symmetry of the model. The bubble in this case is attached to the left boundary, while the location of its right end 
depends on the initial conditions. Finally, on the line p_ = 1~ p+, see Fig.^(b), the bubble coexists with two regions, 
one to the left and one to the right, where the density in both chains is the same. On the left its value is p_, on the 
right it is p+. These regions can coexist, because j{p,p) = j{l — p, 1 — p). In this case, the locations of both ends of 
the bubble depend on the initial conditions. 

In order to verify the correctness of this MF picture, we performed specific Monte-Carlo simulations. In the spirit 
of Derrida et al. (l6| , we introduced two phantom particles A and B designed to track down the left and right end of 
a bubble, respectively. Denoting the position of particle A by a, the rules are: 

+ if ria+i = 0, ma+i = (6) 
a ^ a - I, a {1 - na-i){l ~ nia-i) = (7) 

This means that A moves preferentially to the right in a low-density region and preferentially to the left in a high 
density region, including the bubble region. In contrast to the second-class particles in Ref. [|6|, A does not use the 
sites of the chains and should rather be viewed as moving besides them. Its dynamics therefore does not interfere 
with the that of the chain particles. Analogously, the dynamics of particle B with coordinate b is 

6^6—1, if nf,_i — 1, mb_i — 1 (8) 
6^6+1, if rib-imb-i — (9) 

and it therefore tracks the right end of a bubble. 

One then runs the MC simulations, adds particles A and B at some point and monitors their positions. For 
p_ < 1 — p+ it turns out that A performs a random walk, while B basically sticks to the right boundary, 6 w iV. It 
is then interesting to look at those configurations, where A is at a certain specified site. Thereby one singles out the 
states with a particular size of the bubble region (if there is one) . The average density profile in this case is shown in 
Fig.| for /9_ = 0.15, p+ = 0.77. 
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FIG. 5. Average density profiles as seen from a second-class particle positioned at site 170, provided the density in the first 
chain is higher, for a system with j/=2, p_ = 0.15, p+ = 0.77 and A'' = 400 sites. The average is taken over 6 * 10^ Monte Carlo 
steps after 3 * 10'' steps of equilibration. The horizontal lines indicate the densities p,l — p predicted by equating the currents 
j{p,l — p) = j{p-. , P-), see section tll 
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Quite remarkably, it has the same shape as found in the mean- field calculations |l7|]. This confirms that, indeed, 
there is a very unusual dynamical coexistence of various states. The MC calculations also confirm that the size of the 
bubble region can vary and that the two chains can interchange roles. While this is related to the initial conditions 
in the MF treatment, it happens dynamically here. 

For the case of Fig. ^(c) the situation is different. As mentioned, the system is always half-filled, i.e. pi{t) + p2{t) — 1 
during the evolution. For given overall densities pi, p2, one finds an average profile with one bubble which fills 
essentially the whole system. Thus there is no spatial coexistence of different regions along the chains. However, the 
densities individually fluctuate strongly in an interval 

Pmin < Pia < Pmax (10) 

where roughly pmin = max(/o_, 1 — pmax = 1 — Pmin- A small asymmetry h in the boundary densities (a field in 
magnetic language) suppresses these fluctuations. If one chooses 

P~,i=P--h, p-,2=p-, p+,i=p+, p+,2=P+ + h; (11) 

which makes the chains inequivalent without destroying the particle-hole symmetry, the system locks in at pi — 
Pmin, P2 = Pmax Or vicc vcrsa, depending on the sign of h [Q. Moreover, if one computes formally a 'free energy 
density' f{M) — —log[w{M)]/N at small field, where w{M) is the stationary probability to have M particles in one 
channel, it has exactly the same form as found for a single chain in the vicinity of the line of first-order transitions, 
see Fig. f8 in j^. Therefore, the whole region where this occurs (and which we will determine in more detail below) 
can be viewed as one of first-order transitions. 



IV. STRONG-INTERACTION LIMIT 



The behaviour described at the end of the previous section is found in the whole domain of reservoir densities 

p- < 1/2; p+ > 1/2 (f2) 

if the interaction between the channels is strong, i.e. if ^ f (e <C 1). In this limit, the problem can be discussed 
analytically, as we now show. 

The basic observation is that in this case three of the boundary processes are slow (with rate of order e) and only one 
(where on the left both sites are empty or on the right both sites are full) is fast (rate of order 1). This comes from the 
construction of the boundary rates mentioned above and is discussed in more detail in Appendix A. As a consequence, 
there exist particular configurations, which can be left only with rate e, and which we will call metastable. These are 
such that the system is half filled, with M particles in the first and N ~ M in the second channel, arranged in such a 
way that adjacent sites are always in a particle-hole state. Two examples are shown in Fig. ^. 

o — o — • — o — o — o o — • — o — • — o — • 



FIG. 6. Two examples of metastable states, see text. 



Transitions between such states have the character of avalanches. After a slow process, either at the boundary or 
in the interior, a sequence of N fast processes follows after which a new metastable configuration is reached. If the 
start is in the bulk, four different final configurations with M' = M, M ± I can be reached according to 

with probability 1/2 

with probability 1/4 (13) 
with probability 1/4 

A process in the bulk starts at a wall, where the orientation of the particle-hole pairs changes. For example, the 
first configuration in Fig. |^ has two such walls. For a given Af , the average number of walls can be shown to be 
~ M{1 — M/N). Therefore the average rate for bulk processes from level M is given by 

aM - M(l - M/N)e/4: (14) 
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Similar considerations apply to the processes starting at the boundaries. Here two final configurations M' = M and 
M' = M + 1 or M' = M — 1 can be reached with equal probability. However, the rate does not involve M and thus 
these processes are unimportant compared to the bulk processes except for very small M oy N ~ M . In the following 
we will neglect them. 

The dynamics of the metastable states is thus equivalent to a random walk in the space of occupation numbers 
< M < with a position-dependent rate given by (|l4|). This rate, which in terms of the density pi = p = M/N is 
proportional to p(l — p), is small at the ends of the interval and large in the middle. Therefore the system will spend 
most of its time near M = and M = N . This can be made more precise by writing down the master equation for 
the probability wm to find the value M 

= aM+iWM+i + aM-iWM-i - 2aMWM (15) 

ot 

In a continuum limit, this becomes 

'-^ = $mpMP» (16) 

where D{p) = um/N^ is the diffusion constant. This equation can be related to the associated Legendre differential 
equation. However, for the stationary state it is sufficient to note that the probability current is given by 

/(p) = -^P(pMp)) (17) 

Since this has to be zero, the stationary distribution follows as 

Mp) - p(ilp) ' 0<P<1 (18) 

This simple universal function, which is independent of the parameters of the system, is also seen in simulations 
as illustrated in Fig. 0(a). Actually, the picture shows some deviations at the boundaries of the interval, but there 
the neglected processes should be included. This would lead to finite rates also at p = 0, 1 and prevent w{p) from 
diverging there. In order to normalize it as it stands, one has to leave out a small boundary region of width 6/N. 
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FIG. 7. (a) Stationary probability distribution w{p) for small systems in the seesaw region. Analytical result (11) (continuous 
curve) in comparison with data from Monte-Carlo simulations of systems with v — 8,p- — 0.1, p+ = 0.9 and A'^ = 20,30,50 
sites (crosses, stars and boxes). The systems have evolved for at least 10^ Monte-Carlo steps, the averages are taken over 40 
histories, (b) The corresponding average passage times. 
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Although the system is mainly in configurations with the density in one chain small and in the other large, it is not 
locked into them. The average time in which the two chains interchange roles can be calculated from the first passage 
time formula given in In the continuum case it reads 

l-l-S/N .1/2 , , 

T = / dpw{p) / P = Nln{N/6)/e (19) 

Js/N Js/N D{p')w{p') 

Thus T is smaller than for a homogeneous diffusion process where it varies as N'^. This is due to the large hopping 
rates proportional to N for intermediate values of p. The law is also found numerically, as seen from Fig. Mh). 
However, this result and the considerations so far do not apply to arbitrarily large system sizes. 

If N becomes too large, the typical time 1/Ne after which a new slow process starts, becomes smaller than the time 
T = N needed to complete the first one. Therefore the previous considerations are limited to sizes such that N'^e ^ 1. 
For larger systems, it is not possible to separate the fast and the slow processes. Nevertheless, the system is still 
basically half filled and each time a particle comes in and another particle exits, the number M changes according to 
( p^ ) . Since these processes also determine the flux through the system, their effective rate can be related to j (p, 1 — p) . 
Thus instead of (Q) one has 

aM=j{p,l-p)/2 (20) 

where for j one can use the ring result (^. This should give the behaviour in the thermodynamic limit. Following 
the same steps as before, one then obtains the distribution 

w{p) ~ — -; p within the limits (10) (21) 

j(pA- P) 

for the density p. This is again independent of the boundary values p± , and the interaction enters only via the formula 
for the current. The function is qualitatively similar to ( p^ ) but, according to the form of j (cf. Fig.|^), has a flatter 
shape. 
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FIG. 8. Stationary probability distribution w{p) for a large system in the seesaw phase. Analytical result (14) together with 
with simulation data for a system of A = 300 sites from 1.5- 10^ Monte-Carlo steps, averaging over 40 histories. The parameters 
are f = 4, p- = 0.2 and p+ = 0.8. The result (11) is also shown for comparison (dotted). 



This behaviour is clearly seen in simulations. Fig.g shows numerical results for a system of 300 sites together with 
the analytical prediction. Again, there are deviations in the boundary regions near p = and p — 1 due to the Eq.[l0| 
but apart from that the agreement is very good. The curve (|l|), plotted for comparison, clearly does not fit the data. 
Note that the rate (po| ) does not contain the size any more. Therefore, the passage time T is now proportional to N'^. 
In contrast to the diffusion model treated in P,EO[ there is no exponential increase with the size. 
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V. PHASE DIAGRAM 



The boundaries of the seesaw region discussed above, but also of the mixed region, depend on the interaction 
parameter v and can be found from the following argument. Let us fix p+ > 1/2 and gradually increase p_. For very 
small p_, it takes the time t k, \/p- to fill an empty site from the reservoir, which is longer than the time 1/e for 
slow processes to happen in the bulk. In this case the system will go into the symmetric low-density phase which also 
exists in the lower left corner of the p- — p_|_ phase diagram. As p- < it will tend to minimize its flux according 
to the minimization principle Q and the stationary current will be But this principle, extended to the 

two-channel case, suggests that the low- density phase will become unstable if there is a region (pi, P2) such that the 
currents in both channels are smaller than j (p_ , P-) '■ 



ia{pi,P2) < j{p^,p-). p- < Pl,P2 < P+ 

A transition will take place as soon as the first such point appears. For our system this happens at pi = 1 
and therefore the instability is expected when 

j{p-,p-) = .i{p+A - P+) 



(22) 

P2= P+ 

(23) 



i.e. when the upper and the lower curves in Fig.^ lead to the same current. Actually, according to the results in 
section III, the low-density phase does not vanish completely but coexists with unsymmetric configurations beyond 
that point. 

The part of the region (^2|) along the diagonal pi + p2 = 1 consists of two separated segments (compare Figs.^,c). 
As p^ grows, these segments also grow until they finally merge at a value p* such that 



jip*,p*)=j{l/2,l/2) 



(24) 



From that point onwards, the low-density symmetric region (the presence of which is signalled by the 'cloudy' part 
in Figs.^,c) disappears completely and the seesaw phase described in the previous section takes over. A further 
increase of p_, however, decreases the range within which pi and p2, according to (22), (|o|), may fluctuate. At 
p_ = 1/2 the system crosses the boundary into the symmetric phase with pi = P2 = p+ via a second-order transition 
(see Figs.p^,c). 

These simple arguments are well supported by Monte-Carlo calculations. In Fig. |^ the stationary density distribu- 
tion w(pi, P2) is shown for fixed p+ and four different values of p_. The change of its shape as one crosses the points 
defined by (p|,|3) is clearly seen. 







FIG. 9. Stationary distribution w{p\,p2) as obtained from numerical simulations for a system of 300 sites with = 4, 
p+ = 0.8 and p_ = 0.02, 0.055, 0.06, 0.09 (a-d). The figures show the location of non-zero w, the dark regions within white ones 
correspond to the highest values of w. 
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In Fig. ^ (a) the system is in the conventional low-density phase. The distribution has a boomerang shape due to 
the finite size (incidentally similar to the one found in |^) with the weight however concentrated on a symmetric line 
Pi = P2 near the origin. 

In Fig. ^ (d), on the other hand, there is only weight along the diagonal pi + /32 = 1 which is characteristic of 
the seesaw region. Figures ^ (b) and (c) show the transition region between (23) and ( p^ ) and indicate the mixed 
phase. The 'cloudy' part in Fig.|| (b) and (c) is due to the coexistence of symmetric and asymmetric states mentioned 
above and the increasing weight along the diagonal comes from an increase of density values satisfying Eq. (^) . The 
mixed region, as seen in Figs. |l^(a)-(c), exists for all values e < 1. By contrast, the seesaw region only appears when 
e becomes smaller than tcrit = 1/4, the value where the symmetric current j{p,p) starts to develop a double- hump 
structure. 

The discussion so far has assumed that the value of p^, for which one varies /?_, is not too large. If > p*, where 
p* is defined by (p^, and one increases p_, one will still cross the transition line at the point given by ([2^). With 
further increase of p-, however, the condition ( p^ will not be satisfied, and the system will remain in the mixed 
region. Upon crossing the line p+ = 1 — p-, the coexistence of symmetric it low-density regions with asymmetric 
ones changes into a coexistence of symmetric high-density regions with asymmetric ones. Finally, one ends up in the 
conventional high-density phase. 

One should mention that, if one analyzes the situation more closely with second-class particles, one finds that the 
fraction of symmetric and unsymmetric configurations changes as one moves through the mixed-phase region. This 
leads to the changing 'cloud' mentioned above. However, we do not discuss this in more detail here. 
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FIG. 10. Phase diagram of the model for three different values of the coupling : (a)i^ = 1, (6)i/ — 2, {cjv = 4. Thick (thin) 
lines indicate first (second) order transitions, pi^ax a-nd Pmai denote the positions of the left and right maximum of the curve 
j{p,p) in Fig. 2. 



The complete p_ — p+ -phase diagram is drawn in Fig. 10 (a-c) for three different values of the interaction. The 
sequence shows in particular how the two new regions develop in the upper left corner as v is increased. In the 
non-interacting case, the mixed region shrinks to the upper left part of the diagonal and becomes the single-chain 
transition line. The remaining regions are occupied by symmetric phases with pi — p2, and the boundaries are 
determined by the extremal principle for the current j(p,p), given by (^. This happens because outside the region 
( |l^ ) the evolution of the system at the boundaries is governed by the fast processes. For instance, if p_ > 1/2, the 
injection rate will be fast. This will produce a considerable number of adjacent pairs of particles in the bulk and 
consequently the extraction will also be due to the fast processes. For the symmetric boundary densities which we 
consider, a symmetric bulk situation is to be expected. But then each chain behaves essentially like an independent 
one with current j(p) = j{p,p). Consequently, the problem falls into the class considered in 0. In drawing the 
boundaries, the symmetry of the problem (see the end of section II) was used. They were also checked by simulations. 



VI. CONCLUSION 



We have studied the problem of a two-channel system, coupled to reservoirs of prescribed densities, through which 
a current flows. The aim was to see, which phases one can expect in such a system, and what the principles are which 
govern the transitions between them. The example we took, was a one-parameter model which has a simple stationary 
state on a ring. We found that already this system shows complex behaviour in certain parts of the parameter space. 
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The mixed-phase region is probably the simpler phenomenon. A relatively close analogy in equlibrium statistical 
physics can be found in a system of two ferromagnetic planes which are coupled together antiferromagnetically. The 
coexistence line iJ = 0, T < Tc of the single layers in a uniform field H then widens into a whole region in the T ~ H 
plane, and if one creates domain walls by identical boundary conditions on both layers, these tend to separate in space, 
creating bubble regions of opposite magnetizations, while a non-zero field favours regions of equal magnetizations. 
Thus one finds features as in Figs. 4 and 5. 

The seesaw region with its 'weak' symmetry breaking is more interesting, and we studied it in more detail. One might 
view the phenomenon as a kind of phase- separation between the channels, as opposed to the one along the channels 
in the mixed-phase region. There are also certain similarities to critical phases. On the one hand, the probability 
distribution w{p) resembles that for the order parameter of a finite Ising system at the critical temperature One 
the other hand, its properties do not depend (except for the left and right limits) on the boundary values /o±, a feature 
which it shares with the critical (maximal-current) phase found in the present model. 

One should note that our model differs from the lattice-gas models studied in since the particles cannot hop 

between the chains. It would be interesting to see whether relaxing these conservation laws changes the situation 
strongly. Leaving this aspect aside, however, it seems that the results are rather general. We did use the microscopic 
details for arguing in the paper, but the final results do not depend on them directly. It is only the fiux, which 
determines the phase transition lines and characterizes the new phases. Thus one can hope that pursuing this approach 
would allow to formulate generic principles which govern the phase transitions in such multi-channel systems. 
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APPENDIX A: BOUNDARY RATES 

In our model we view the reservoirs as extensions of the system having the same properties as the ring . This 
permits a natural definition of the boundary rates. Consider, for example, the first process in Fig. |l|. On the ring, it 
leads to the average current 

{ji) = a(n/(l - m/))(n;+i(l - ?7i/+i)) (Al) 

between sites I and I + 1. Here the product form of the stationary state has been used. If site I belongs to the left 
reservoir, one combines a with the first factor, evaluated for the prescribed boundary densities pi ~ p-i ~ This 
gives the rate — a{n{l — m))p-. Analogously, if site (I + 1) belongs to the right reservoir, one combines a with 
the second factor and defines a+ — a({l — n){l — m))p^. 

In this way, one obtains boundary rates which are determined by the bulk ones, multiplied by the weight of the 
boundary configuration involved in the process. With a = (3 — "f — I, one has 

a^=P--P-; P-=l-=P-\ e-=e(P--P-); (A2) 
a+ = 7+ =P+ + 1 - 2p+; 13+ = p+ - p+; e+ = e{p+ ~ p+) (A3) 

where p± = p{p±) = {nm)p± and p{p) is given by 

P{p) ^ P - 



{2p-lf -Aep{l~p)-1 



/2{l-e) (A4) 



For e ^ 1, configurations with two particles (holes) at one site have a small weight if p < 1/2 (p > 1/2). Thus, for 
/9_ < 1/2 and p+ > 1/2, the rates 7_, a+, 7+ are proportional to e, as are the rates e±. In this case, only one fast 
boundary process remains at each end. 

With these quantities one can then write down the total rates for the boundary processes. For example, if the first 
sites of both chains are empty, the rate with which a particle enters one chain is given by («_ +7-)- As a consequence, 
the currents at the left and the right end of the first chain are, respectively 

jl = (a_ + 7_) (1 - 7ii) (1 - mi) + (/?_ + e_) (1 - ni) mi; (A5) 
j]_ = (/3+ + 7+) nNiriN + (a+ + e+) (1 - mAr) ; (A6) 
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If p- = p+, the stationary state of the system is automatically the same as for the ring, i.e. the density is constant 
everywhere. This is the basic motivation for the approach described here. For e = 1 everything reduces to the 
single-chain problem [||J|] . 



APPENDIX B: MEAN-FIELD EQUATIONS 



In the mean-field approximation, one neglects the correlations between different sites. In our case, the system on 
a ring has no correlations in the stationary state {nink) — {ni){nk). Two adjacent sites on different chains, however, 
are correlated. Therefore one should take the product pi = nimi as an independent variable. Keeping this in mind, 
the mean-field equations resulting from the gain and the loss processes become (for a = /3 = 7 = 1 and leaving out 
the averaging brackets for simplicity) 



duk 
dt 

dpk 
dt 



— Jk Jk+1 

[nfe_i(l - Uk) + (e - l){nk-i ~ Pk-i)(mk ~ Pk)] ~ ["-fe(l - '^fc+i) + (e - l)("-fe - Pk){mk+i ~ Pk+i)] (Bl) 
{nk -Pk)iemk-i + (1 - e)Pfe-i) + ("T-fe - Pk){enk-i + (1 - e)Pfe-i) -Pfe(2 - rik+i - mk+i). (B2) 



The equation for the rrifc is obtained by substituting m ^ n in (Bl). 

The homogeneous solution of these equations = rn, Uk = n,pk = p leads to a quadratic equation for p. When 
substituted into the expression for the current in (Bl), this reproduces the value (^ because the mean- field equations 
are exact in this case. 

For the open system, the equations (Bl B2) should be supplemented by the boundary conditions tiq = niQ = 
/3_,n7v+i — m^v+i — p+. For pk one takes the homogeneous solution at both ends. 
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